##############
## American Identity, Status Threat, and Backlash - Study 3, Public
##############

##############
## Load data
##############
rm(list=ls())
S3_D1 <- read.csv("Study 3 raw.csv", header = TRUE)
head(S3_D1)
nrow(S3_D1)
ncol(S3_D1)

##############
## Exclude participants who failed the attention check
##############
variable.names(S3_D1)
table(S3_D1$Att_Check2)
nrow(S3_D1)

## remove based on main attention check
S3_D2 <- S3_D1[which((S3_D1$Att_Check2) %in% c("meow","Meow","meow meow","meow =^.^=","meow!","meow or purr","Mhew :)","Mew","purr","No noise cat made","meows","meow!","Meow.","no","mew or meow, depending on how silly the cat is.",
                                                    "Meow man","pur","MEOW","meow or miaow","meeow","meeowww","purrrrrrr","Meow meow","Meow!","moi","meows","meyow","none","MEOWW","Miya miya","meowwww","MEOW!","meow purrrrrrrrr",
                                                    "meow!","mew","mewow","miau"," Meow!","this cat looks irritated, so it would be a growly kind of meow. Give him a treat.",
                                                    "PURR", "purrrrr", "Purr","Mrrp?  Mrrrrrow.  mrrrrAAAAAA!","miow","miaow", " Miao or meow","Meow  Prrr Prrr","meow (or purrrrrr...) ","meow (or nya in japanese) ", "meo", "mea", "Meows",
                                                    "MEOOOOOW", "Bark", "meow (or nya in japanese)","meow (or purrrrrr...)","Either a purr or a meow", "Meaw", "meow, hiss!","meow, meow, meeeeooww, mew, meow","meow and purr", "Meow Hiss",
                                                    "meo","meoww","myaw","meeyow","Miao or meow","Meawww","meowwwwww","Meoww","meow/purr", "cat", "meow", "MEOW MEOW", "Meow Meow", " meow", "meoow",
                                                    "MEOW MEOW", "Meow Meow","purr purr meow hissss","meow ","Mew, meow...","MEOW, purrrrrr!","	meooow!","meooow!","Meow","Meow, hiss","myau","Purr! :)","Meowww","Meow ","MEOOOWW","Meowwwew","Meow! ","Meowwww","Meow meow! =^.^=","Meow ",
                                               "meooooooooooooooow","Purrr","Meow. Or maybe that prrt noise cats make when they're saying hello.","meow.","meow/hiss","meow purr?","mey","Meow! Purr purr! Hiss!","Meow miau","Purring or meow","meow prrr","purrrr","meooooooooooooooow","Purrr",
                                               "miaw")), ]
nrow(S3_D2)

## check those who failed
S3_D0 <- S3_D1[-which((S3_D1$Att_Check2) %in% c("meow","Meow","meow meow","meow =^.^=","meow!","meow or purr","Mhew :)","Mew","purr","No noise cat made","meows","meow!","Meow.","no","mew or meow, depending on how silly the cat is.",
                                                "Meow man","pur","MEOW","meow or miaow","meeow","meeowww","purrrrrrr","Meow meow","Meow!","moi","meows","meyow","none","MEOWW","Miya miya","meowwww","MEOW!","meow purrrrrrrrr",
                                                "meow!","mew","mewow","miau"," Meow!","this cat looks irritated, so it would be a growly kind of meow. Give him a treat.",
                                                "PURR", "purrrrr", "Purr","Mrrp?  Mrrrrrow.  mrrrrAAAAAA!","miow","miaow", " Miao or meow","Meow  Prrr Prrr","meow (or purrrrrr...) ","meow (or nya in japanese) ", "meo", "mea", "Meows",
                                                "MEOOOOOW", "Bark", "meow (or nya in japanese)","meow (or purrrrrr...)","Either a purr or a meow", "Meaw", "meow, hiss!","meow, meow, meeeeooww, mew, meow","meow and purr", "Meow Hiss",
                                                "meo","meoww","myaw","meeyow","Miao or meow","Meawww","meowwwwww","Meoww","meow/purr", "cat", "meow", "MEOW MEOW", "Meow Meow", " meow", "meoow",
                                                "MEOW MEOW", "Meow Meow","purr purr meow hissss","meow ","Mew, meow...","MEOW, purrrrrr!","	meooow!","meooow!","Meow","Meow, hiss","myau","Purr! :)","Meowww","Meow ","MEOOOWW","Meowwwew","Meow! ","Meowwww","Meow meow! =^.^=","Meow ",
                                                "meooooooooooooooow","Purrr","Meow. Or maybe that prrt noise cats make when they're saying hello.","meow.","meow/hiss","meow purr?","mey","Meow! Purr purr! Hiss!","Meow miau","Purring or meow","meow prrr","purrrr","meooooooooooooooow","Purrr",
                                                "miaw")), ]
nrow(S3_D0)
#view(S3_D0$Att_Check2)
library(rstatix)
chisq_test(S3_D0$Att_Check2, S3_D0$Condition) ## inattentiveness does not vary by condition


## remove based on second attention check
table(S3_D2$White_ID_6)
S3_D3 <- subset(S3_D2, S3_D2$White_ID_6=="3") 
nrow(S3_D3)
chisq_test(S3_D2$White_ID_6, S3_D2$Condition) ## inattentiveness does not vary by condition


## Exclude all non-whites
table(S3_D3$Ethnicity)
S3_D3 <- subset(S3_D3, S3_D3$Ethnicity==1) 
nrow(S3_D3)

##############
## Demographic summary + balance tests
##############
variable.names(S3_D3)

S3_D3$Condition [S3_D3$Condition == 1] <- "Control"
S3_D3$Condition [S3_D3$Condition == 2] <- "Common Identity"

is.factor(S3_D3$Condition)

## Demographics
mean(S3_D3$Age, na.rm = TRUE)
sd(S3_D3$Age, na.rm = TRUE)
round(table(S3_D3$Ideology)/nrow(S3_D3)*100,2) ## percentage liberal
round(table(S3_D3$Gender)/nrow(S3_D3)*100,2) ## percentage female

## Balance tests
nrow(S3_D3)
table(S3_D3$Condition)
chisq.test(table(S3_D3$Condition, S3_D3$Age))
chisq.test(table(S3_D3$Condition, S3_D3$Gender))
chisq.test(table(S3_D3$Condition, S3_D3$Education))
chisq.test(table(S3_D3$Condition, S3_D3$Ideology))


##############
## Recode reverse-coded items and rename variables
##############
## recode second status threat item
library(car)
S3_D3$Status_Threat_2 = car::recode(S3_D3$Status_Threat_2, "1=7; 2=6; 3=5; 4=4; 5=3; 6=2; 7=1; else = NA")
S3_D3$Status_Threat_4 = car::recode(S3_D3$Status_Threat_4, "1=7; 2=6; 3=5; 4=4; 5=3; 6=2; 7=1; else = NA")

table(S3_D3$American_ID_4)
table(S3_D3$Americanness_1)
S3_D3$American_ID_1 = car::recode(S3_D3$American_ID_1, "32=1; 33=2; 34=3; 35=4;36=5; 37=6; 38=7; else = NA")
S3_D3$American_ID_2 = car::recode(S3_D3$American_ID_2, "32=1; 33=2; 34=3; 35=4;36=5; 37=6; 38=7; else = NA")
S3_D3$American_ID_3 = car::recode(S3_D3$American_ID_3, "32=1; 33=2; 34=3; 35=4;36=5; 37=6; 38=7; else = NA")
S3_D3$American_ID_4 = car::recode(S3_D3$American_ID_4, "32=1; 33=2; 34=3; 35=4;36=5; 37=6; 38=7; else = NA")

S3_D3$Americanness_1 = car::recode(S3_D3$Americanness_1, "32=1; 33=2; 34=3; 35=4;36=5; 37=6; 38=7; else = NA")
S3_D3$Americanness_2 = car::recode(S3_D3$Americanness_2, "32=1; 33=2; 34=3; 35=4;36=5; 37=6; 38=7; else = NA")
S3_D3$Americanness_3 = car::recode(S3_D3$Americanness_3, "32=1; 33=2; 34=3; 35=4;36=5; 37=6; 38=7; else = NA")
S3_D3$Americanness_4 = car::recode(S3_D3$Americanness_4, "32=1; 33=2; 34=3; 35=4;36=5; 37=6; 38=7; else = NA")
S3_D3$Americanness_5 = car::recode(S3_D3$Americanness_5, "32=1; 33=2; 34=3; 35=4;36=5; 37=6; 38=7; else = NA")
S3_D3$Americanness_6 = car::recode(S3_D3$Americanness_6, "32=1; 33=2; 34=3; 35=4;36=5; 37=6; 38=7; else = NA")

##############
## Factorial structure outgroup thermometers
##############
S3_D3_thermometers <- S3_D3[,c(32:36)]
head(S3_D3_thermometers)

## Visualize with scree plot
library(ggplot2)
library(psych)
fafitfree <- fa(S3_D3_thermometers,nfactors = 1, rotate = "promax", fm = "ml")
n_factors <- length(fafitfree$e.values)
scree     <- data.frame(
  Factor_n =  as.factor(1:n_factors), 
  Eigenvalue = fafitfree$e.values)
ggplot(scree, aes(x = Factor_n, y = Eigenvalue, group = 1)) + 
  geom_point() + geom_line() +
  xlab("Number of factors") +
  ylab("Initial eigenvalue") +
  labs( title = "Scree Plot", 
        subtitle = "(Based on the unreduced correlation matrix)")


##############
## Factorial structure outgroup Americanness
##############
S3_D3_americanness <- S3_D3[,c(40:42)]
head(S3_D3_americanness)

## Visualize with scree plot
library(ggplot2)
fafitfree <- fa(S3_D3_americanness,nfactors = 1, rotate = "promax", fm = "ml")
n_factors <- length(fafitfree$e.values)
scree     <- data.frame(
  Factor_n =  as.factor(1:n_factors), 
  Eigenvalue = fafitfree$e.values)
ggplot(scree, aes(x = Factor_n, y = Eigenvalue, group = 1)) + 
  geom_point() + geom_line() +
  xlab("Number of factors") +
  ylab("Initial eigenvalue") +
  labs( title = "Scree Plot", 
        subtitle = "(Based on the unreduced correlation matrix)")


##############
## Factorial structure National Pride
##############
S3_D3_pride <- S3_D3[,c(22:28)]
head(S3_D3_pride)

## Visualize with scree plot
library(ggplot2)
fafitfree <- fa(S3_D3_pride,nfactors = 1, rotate = "promax", fm = "ml")
n_factors <- length(fafitfree$e.values)
scree     <- data.frame(
  Factor_n =  as.factor(1:n_factors), 
  Eigenvalue = fafitfree$e.values)
ggplot(scree, aes(x = Factor_n, y = Eigenvalue, group = 1)) + 
  geom_point() + geom_line() +
  xlab("Number of factors") +
  ylab("Initial eigenvalue") +
  labs( title = "Scree Plot", 
        subtitle = "(Based on the unreduced correlation matrix)")


##############
## EFA Conceptions of nationhood
##############
S3_D3_nationhood <- S3_D3[,c(14:28,46:50)]
head(S3_D3_nationhood)

## Visualize with scree plot
library(ggplot2)
fafitfree <- fa(S3_D3_nationhood,nfactors = 1, rotate = "promax", fm = "ml")
n_factors <- length(fafitfree$e.values)
scree     <- data.frame(
  Factor_n =  as.factor(1:n_factors), 
  Eigenvalue = fafitfree$e.values)
ggplot(scree, aes(x = Factor_n, y = Eigenvalue, group = 1)) + 
  geom_point() + geom_line() +
  xlab("Number of factors") +
  ylab("Initial eigenvalue") +
  labs( title = "Scree Plot", 
        subtitle = "(Based on the unreduced correlation matrix)")


## KMO-test
KMO(S3_D3_nationhood) ## overall MSA = .93, so high to 1. Data thus suitable for FA

## Bartlett test
cortest.bartlett(S3_D3_nationhood) ## showing variables are correlated, H0 rejected. Data is not an identity matrix and FA is thus supported.

## scree-test
scree(S3_D3_nationhood, pc=FALSE) ## suggesting several factors

## parallel analysis
fa.parallel(S3_D3_nationhood, fa="fa") ## suggesting several factors

## EFA (factanal)
S3_D3_nationhood2 <- S3_D3_nationhood[complete.cases(S3_D3_nationhood), ] ## choose only complete cases
nrow(S3_D3_nationhood2)
S3_D3_nationhood3 <- factanal(S3_D3_nationhood2, 5, rotation="promax")
print(S3_D3_nationhood3, cutof = .300)

load_S3_D3_nationhood3 <- S3_D3_nationhood3$loadings[,1:5] ## plot S3_D3_nationhood3
plot(load_S3_D3_nationhood3,type="n") # set up plot
text(load_S3_D3_nationhood3,labels=names(S3_D3_nationhood2),cex=.7)

loads_S3_D3_nationhood3 <- S3_D3_nationhood3$loadings ## visualize factor model
fa.diagram(loads_S3_D3_nationhood3)

##############
## CFA Conceptions of nationhood
##############
library(lavaan)
CFA_Conceptions_S3 <- '
## Direct effect condition on outcomes (path c) 
AM_ID =~ American_ID_1 + American_ID_2 + American_ID_3 + American_ID_4
White_ID =~ White_ID_1 + White_ID_2 + White_ID_3 + White_ID_4 + White_ID_5
Patr =~ Patriotism_1 + Patriotism_2
Natio =~ Nationalism_1 + Nationalism_2
Pride =~ Pride_1 + Pride_2 + Pride_3 + Pride_4 + Pride_5 + Pride_6 + Pride_7

AM_ID ~~ White_ID
AM_ID ~~ Patr
AM_ID ~~ Natio
AM_ID ~~ Pride
White_ID ~~ Patr
White_ID ~~ Natio
White_ID ~~ Pride
Patr ~~ Natio
Patr ~~ Pride
Natio ~~ Pride

## modindices
Pride_3 ~~ Pride_6
American_ID_4 ~~ White_ID_5
'

CFA_Conceptions_S3_fit <- cfa(CFA_Conceptions_S3, se = "bootstrap", bootstrap = 100, data = S3_D3)
summary(CFA_Conceptions_S3_fit, fit.measures = TRUE, rsquare = TRUE)
modindices(CFA_Conceptions_S3_fit, max = 30, sort = TRUE)

##############
## Reliabilities
##############
colnames(S3_D3)
which( colnames(S3_D3)=="Pride_7" )

psych::alpha(select(S3_D3, 9:13)) ## Status threat
psych::alpha(select(S3_D3, 14:17)) ## American ID
psych::alpha(select(S3_D3, 46:50)) ## White ID
psych::alpha(select(S3_D3, 18:19)) ## Symbolic patriotism
psych::alpha(select(S3_D3, 20:21)) ## Nationalism
psych::alpha(select(S3_D3, 22:28)) ## National pride

psych::alpha(select(S3_D3, 32:36)) ## Outgroup thermometers
psych::alpha(select(S3_D3, 39:41)) ## Perceived Americanness


##############
## Mean Scales
##############
S3_D3$mean_Status_Threat <- (S3_D3$Status_Threat_1 + S3_D3$Status_Threat_2 + S3_D3$Status_Threat_3 + S3_D3$Status_Threat_4 + S3_D3$Status_Threat_5)/5
S3_D3$mean_American_Identity <- (S3_D3$American_ID_1 + S3_D3$American_ID_2 + S3_D3$American_ID_3 + S3_D3$American_ID_4)/4
S3_D3$mean_White_Identity <- (S3_D3$White_ID_1 + S3_D3$White_ID_2 + S3_D3$White_ID_3 + S3_D3$White_ID_4 + S3_D3$White_ID_5)/5
S3_D3$mean_Symbolic_Patriotism <- (S3_D3$Patriotism_1 + S3_D3$Patriotism_2)/2
S3_D3$mean_Nationalism <- (S3_D3$Nationalism_1 + S3_D3$Nationalism_2)/2
S3_D3$mean_National_Pride <- (S3_D3$Pride_1 + S3_D3$Pride_2 + S3_D3$Pride_3 + S3_D3$Pride_4 + S3_D3$Pride_5 + S3_D3$Pride_6 + S3_D3$Pride_7)/7
S3_D3$mean_Outgroup_Warmth <- (S3_D3$Thermometers_6 + S3_D3$Thermometers_7 + S3_D3$Thermometers_8 + S3_D3$Thermometers_9 + S3_D3$Thermometers_10)/5
S3_D3$mean_Perceived_Americanness <- (S3_D3$Americanness_4 + S3_D3$Americanness_5 + S3_D3$Americanness_6)/3

##############
## Rescaling
##############
S3_D3$mean_Status_Threat_scaled <- (S3_D3$mean_Status_Threat - min(S3_D3$mean_Status_Threat, na.rm = TRUE)) / (max(S3_D3$mean_Status_Threat, na.rm = TRUE) - min(S3_D3$mean_Status_Threat, na.rm = TRUE))
S3_D3$mean_American_Identity_scaled <- (S3_D3$mean_American_Identity - min(S3_D3$mean_American_Identity, na.rm = TRUE)) / (max(S3_D3$mean_American_Identity, na.rm = TRUE) - min(S3_D3$mean_American_Identity, na.rm = TRUE))
S3_D3$mean_White_Identity_scaled <- (S3_D3$mean_White_Identity - min(S3_D3$mean_White_Identity, na.rm = TRUE)) / (max(S3_D3$mean_White_Identity, na.rm = TRUE) - min(S3_D3$mean_White_Identity, na.rm = TRUE))
S3_D3$mean_Symbolic_Patriotism_scaled <- (S3_D3$mean_Symbolic_Patriotism - min(S3_D3$mean_Symbolic_Patriotism, na.rm = TRUE)) / (max(S3_D3$mean_Symbolic_Patriotism, na.rm = TRUE) - min(S3_D3$mean_Symbolic_Patriotism, na.rm = TRUE))
S3_D3$mean_Nationalism_scaled <- (S3_D3$mean_Nationalism - min(S3_D3$mean_Nationalism, na.rm = TRUE)) / (max(S3_D3$mean_Nationalism, na.rm = TRUE) - min(S3_D3$mean_Nationalism, na.rm = TRUE))
S3_D3$mean_National_Pride_scaled <- (S3_D3$mean_National_Pride - min(S3_D3$mean_National_Pride, na.rm = TRUE)) / (max(S3_D3$mean_National_Pride, na.rm = TRUE) - min(S3_D3$mean_National_Pride, na.rm = TRUE))
S3_D3$mean_Outgroup_Warmth_scaled <- (S3_D3$mean_Outgroup_Warmth - min(S3_D3$mean_Outgroup_Warmth, na.rm = TRUE)) / (max(S3_D3$mean_Outgroup_Warmth, na.rm = TRUE) - min(S3_D3$mean_Outgroup_Warmth, na.rm = TRUE))
S3_D3$IOS_scaled <- (S3_D3$IOS - min(S3_D3$IOS, na.rm = TRUE)) / (max(S3_D3$IOS, na.rm = TRUE) - min(S3_D3$IOS, na.rm = TRUE))
S3_D3$mean_Perceived_Americanness_scaled <- (S3_D3$mean_Perceived_Americanness - min(S3_D3$mean_Perceived_Americanness, na.rm = TRUE)) / (max(S3_D3$mean_Perceived_Americanness, na.rm = TRUE) - min(S3_D3$mean_Perceived_Americanness, na.rm = TRUE))
S3_D3$Dems_scaled <- (S3_D3$Thermometers_3 - min(S3_D3$Thermometers_3, na.rm = TRUE)) / (max(S3_D3$Thermometers_3, na.rm = TRUE) - min(S3_D3$Thermometers_3, na.rm = TRUE))
S3_D3$Reps_scaled <- (S3_D3$Thermometers_4 - min(S3_D3$Thermometers_4, na.rm = TRUE)) / (max(S3_D3$Thermometers_4, na.rm = TRUE) - min(S3_D3$Thermometers_4, na.rm = TRUE))

S3_D3$Blacks_scaled <- (S3_D3$Thermometers_6 - min(S3_D3$Thermometers_6, na.rm = TRUE)) / (max(S3_D3$Thermometers_6, na.rm = TRUE) - min(S3_D3$Thermometers_6, na.rm = TRUE))
S3_D3$Asians_scaled <- (S3_D3$Thermometers_7 - min(S3_D3$Thermometers_7, na.rm = TRUE)) / (max(S3_D3$Thermometers_7, na.rm = TRUE) - min(S3_D3$Thermometers_7, na.rm = TRUE))
S3_D3$Hispanics_scaled <- (S3_D3$Thermometers_8 - min(S3_D3$Thermometers_8, na.rm = TRUE)) / (max(S3_D3$Thermometers_8, na.rm = TRUE) - min(S3_D3$Thermometers_8, na.rm = TRUE))
S3_D3$Illegal_scaled <- (S3_D3$Thermometers_9 - min(S3_D3$Thermometers_9, na.rm = TRUE)) / (max(S3_D3$Thermometers_9, na.rm = TRUE) - min(S3_D3$Thermometers_9, na.rm = TRUE))
S3_D3$Legal_scaled <- (S3_D3$Thermometers_10 - min(S3_D3$Thermometers_10, na.rm = TRUE)) / (max(S3_D3$Thermometers_10, na.rm = TRUE) - min(S3_D3$Thermometers_10, na.rm = TRUE))
S3_D3$Americanness_Blacks_scaled <- (S3_D3$Americanness_4 - min(S3_D3$Americanness_4, na.rm = TRUE)) / (max(S3_D3$Americanness_4, na.rm = TRUE) - min(S3_D3$Americanness_4, na.rm = TRUE))
S3_D3$Americanness_Asians_scaled <- (S3_D3$Americanness_5 - min(S3_D3$Americanness_5, na.rm = TRUE)) / (max(S3_D3$Americanness_5, na.rm = TRUE) - min(S3_D3$Americanness_5, na.rm = TRUE))
S3_D3$Americanness_Hispanics_scaled <- (S3_D3$Americanness_6 - min(S3_D3$Americanness_6, na.rm = TRUE)) / (max(S3_D3$Americanness_6, na.rm = TRUE) - min(S3_D3$Americanness_6, na.rm = TRUE))
S3_D3$Americanness_Whites_scaled <- (S3_D3$Americanness_3 - min(S3_D3$Americanness_3, na.rm = TRUE)) / (max(S3_D3$Americanness_3, na.rm = TRUE) - min(S3_D3$Americanness_3, na.rm = TRUE))
S3_D3$Americanness_Dems_scaled <- (S3_D3$Americanness_1 - min(S3_D3$Americanness_1, na.rm = TRUE)) / (max(S3_D3$Americanness_1, na.rm = TRUE) - min(S3_D3$Americanness_1, na.rm = TRUE))
S3_D3$Americanness_Reps_scaled <- (S3_D3$Americanness_2 - min(S3_D3$Americanness_2, na.rm = TRUE)) / (max(S3_D3$Americanness_2, na.rm = TRUE) - min(S3_D3$Americanness_2, na.rm = TRUE))

S3_D3$Ideology_scaled <- (S3_D3$Ideology - min(S3_D3$Ideology, na.rm = TRUE)) / (max(S3_D3$Ideology, na.rm = TRUE) - min(S3_D3$Ideology, na.rm = TRUE))

##############
## Correlations Conceptions of nationhood
##############
library(corrplot)
library(dplyr)
library(apaTables)

correlations_constructs_S3 <- as.data.frame(cbind(S3_D3$mean_American_Identity_scaled, S3_D3$mean_White_Identity_scaled,S3_D3$mean_Symbolic_Patriotism_scaled,S3_D3$mean_Nationalism_scaled,S3_D3$mean_National_Pride_scaled,
                                                  S3_D3$mean_Status_Threat_scaled))
correlations_constructs_S3_2 <- correlations_constructs_S3 %>% 
  rename(
    "American ID" = V1,
    "White ID" = V2,    
    "Symbolic patriotism" = V3,
    "Nationalism" = V4,
    "National pride" = V5,
    "Status Threat" = V6)
str(correlations_constructs_S3_2)

## Correlation Table
apa.cor.table(correlations_constructs_S3_2, table.number = 2, show.sig.stars = TRUE, show.conf.interval = FALSE, filename = "construct_correlations.doc")
cortable_S3 <- apa.cor.table(correlations_constructs_S3_2, table.number = 2, show.sig.stars = TRUE, show.conf.interval = FALSE, filename = "construct_correlations.doc")

## Significance test of correlations
correlations_constructs_S3_sig <- cor.mtest(correlations_constructs_S3_2)

## Correlogram 
construct_corr_S3 <- cor(correlations_constructs_S3_2, use = "complete.obs")
head(round(construct_corr_S3,2))

corrplot(construct_corr_S3, method="color", addCoef.col = "black", p.mat = correlations_constructs_S3_sig$p, sig.level = 0.001, type = "lower", insig = "blank", diag = TRUE,
         title = "",
         tl.cex = 0.5, number.cex = 0.4, cl.align.text = "l", tl.pos = "ld", tl.col="black", tl.offset = 3)
ggsave("Correlations_S3.jpeg", scale = 1, dpi="retina", dev='jpeg', height=5, width=5, units="in")


##############
## Descriptives
##############
table(S3_D3$Condition)
S3_D3$Condition <- factor(S3_D3$Condition, levels = c("Control", "Common Identity"))
S3_D3$Condition_num <- as.numeric(as.factor(S3_D3$Condition))
table(S3_D3$Condition_num)

S3_t1 <- t(round(S3_D3 %>%
                   group_by(Condition_num) %>%
                   summarise_at(vars(c(mean_Status_Threat_scaled, 
                                       mean_American_Identity_scaled, mean_White_Identity_scaled,mean_Symbolic_Patriotism_scaled,mean_Nationalism_scaled,mean_National_Pride_scaled,
                                       mean_Outgroup_Warmth_scaled, IOS_scaled,mean_Perceived_Americanness_scaled,
                                       Dems_scaled,Reps_scaled,
                                       Blacks_scaled,Asians_scaled,Hispanics_scaled,Illegal_scaled,Legal_scaled,
                                       Americanness_Blacks_scaled,Americanness_Asians_scaled,Americanness_Hispanics_scaled,
                                       Americanness_Whites_scaled,Americanness_Dems_scaled,Americanness_Reps_scaled)), list(mean = mean, sd = sd), na.rm = TRUE),2))
write.table(S3_t1, file = "S3_t1.txt", sep = ",", quote = FALSE, row.names = F)


##############
## Simple mean differences (descriptive, not used for hypothesis testing)
##############
## Status threat
S3_aov1 <- aov(mean_Status_Threat_scaled ~ factor(Condition), data = S3_D3)
summary(S3_aov1)
S3_reg1 <- lm(mean_Status_Threat_scaled ~ factor(Condition), data = S3_D3)
summary(S3_reg1)

library(effectsize)
cohens_f(S3_reg1, partial = FALSE)

## American identity
S3_aov2 <- aov(mean_American_Identity_scaled ~ factor(Condition), data = S3_D3)
summary(S3_aov2)
S3_reg2 <- lm(mean_American_Identity_scaled ~ factor(Condition), data = S3_D3)
summary(S3_reg2)
cohens_f(S3_reg2, partial = FALSE)

## White identity
S3_aov3 <- aov(mean_White_Identity_scaled ~ factor(Condition), data = S3_D3)
summary(S3_aov3)
S3_reg3 <- lm(mean_White_Identity_scaled ~ factor(Condition), data = S3_D3)
summary(S3_reg3)
cohens_f(S3_reg3, partial = FALSE)


## Symbolic patriotism
S3_aov4 <- aov(mean_Symbolic_Patriotism_scaled ~ factor(Condition), data = S3_D3)
summary(S3_aov4)
S3_reg4 <- lm(mean_Symbolic_Patriotism_scaled ~ factor(Condition), data = S3_D3)
summary(S3_reg4)
cohens_f(S3_reg4, partial = FALSE)

## Nationalism
S3_aov5 <- aov(mean_Nationalism_scaled ~ factor(Condition), data = S3_D3)
summary(S3_aov5)
S3_reg5 <- lm(mean_Nationalism_scaled ~ factor(Condition), data = S3_D3)
summary(S3_reg5)
cohens_f(S3_reg5, partial = FALSE)


## National pride
S3_aov6 <- aov(mean_National_Pride_scaled ~ factor(Condition), data = S3_D3)
summary(S3_aov6)
S3_reg6 <- lm(mean_National_Pride_scaled ~ factor(Condition), data = S3_D3)
summary(S3_reg6)
cohens_f(S3_reg6, partial = FALSE)


## Outgroup warmth
S3_aov7 <- aov(mean_Outgroup_Warmth_scaled ~ factor(Condition), data = S3_D3)
summary(S3_aov7)
S3_reg7 <- lm(mean_Outgroup_Warmth_scaled ~ factor(Condition), data = S3_D3)
summary(S3_reg7)

## Inclusiveness
S3_aov8 <- aov(IOS_scaled ~ factor(Condition), data = S3_D3)
summary(S3_aov8)
S3_reg8 <- lm(IOS_scaled ~ factor(Condition), data = S3_D3)
summary(S3_reg8)

## Americanness
S3_aov9 <- aov(mean_Perceived_Americanness_scaled ~ factor(Condition), data = S3_D3)
summary(S3_aov9)
S3_reg9 <- lm(mean_Perceived_Americanness_scaled ~ factor(Condition), data = S3_D3)
summary(S3_reg9)

## Democrat warmth
S3_aov10 <- aov(Dems_scaled ~ factor(Condition), data = S3_D3)
summary(S3_aov10)
S3_reg10 <- lm(Dems_scaled ~ factor(Condition), data = S3_D3)
summary(S3_reg10)

## Republican warmth
S3_aov11 <- aov(Reps_scaled ~ factor(Condition), data = S3_D3)
summary(S3_aov11)
S3_reg11 <- lm(Reps_scaled ~ factor(Condition), data = S3_D3)
summary(S3_reg11)

## Separate outgroups: thermometers
S3_aov12a <- aov(Blacks_scaled ~ factor(Condition), data = S3_D3)
S3_aov12b <- aov(Asians_scaled ~ factor(Condition), data = S3_D3)
S3_aov12c <- aov(Hispanics_scaled ~ factor(Condition), data = S3_D3)
S3_aov12d <- aov(Illegal_scaled ~ factor(Condition), data = S3_D3)
S3_aov12e <- aov(Legal_scaled ~ factor(Condition), data = S3_D3)
summary(S3_aov12a)
summary(S3_aov12b)
summary(S3_aov12c)
summary(S3_aov12d)
summary(S3_aov12e)

## Separate outgroups: Americanness
S3_aov13a <- aov(Americanness_Blacks_scaled ~ factor(Condition), data = S3_D3)
S3_aov13b <- aov(Americanness_Asians_scaled ~ factor(Condition), data = S3_D3)
S3_aov13c <- aov(Americanness_Hispanics_scaled ~ factor(Condition), data = S3_D3)
S3_aov13d <- aov(Americanness_Whites_scaled ~ factor(Condition), data = S3_D3)
summary(S3_aov13a)
summary(S3_aov13b)
summary(S3_aov13c)
summary(S3_aov13d)


##############
## Moderating ideology
##############
S3_aov1_ideo <- aov(mean_Status_Threat_scaled ~ factor(Condition)*Ideology_scaled, data = S3_D3)
summary(S3_aov1_ideo)
S3_reg1_ideo <- lm(mean_Status_Threat_scaled ~ factor(Condition)*Ideology_scaled, data = S3_D3)
summary(S3_reg1_ideo)

S3_aov2_ideo <- aov(mean_American_Identity_scaled ~ factor(Condition)*Ideology_scaled, data = S3_D3)
summary(S3_aov2_ideo)
S3_reg2_ideo <- lm(mean_American_Identity_scaled ~ factor(Condition)*Ideology_scaled, data = S3_D3)
summary(S3_reg2_ideo)

S3_aov3_ideo <- aov(mean_Status_Threat_scaled ~ mean_American_Identity_scaled*Ideology_scaled, data = S3_D3)
summary(S3_aov3_ideo)
S3_reg3_ideo <- lm(mean_Status_Threat_scaled ~ mean_American_Identity_scaled*Ideology_scaled, data = S3_D3)
summary(S3_reg3_ideo)

S3_reg4_ideo <- lm(mean_White_Identity_scaled ~ factor(Condition)*Ideology_scaled, data = S3_D3)
S3_reg5_ideo <- lm(mean_Symbolic_Patriotism_scaled ~ factor(Condition)*Ideology_scaled, data = S3_D3)
S3_reg6_ideo <- lm(mean_Nationalism_scaled ~ factor(Condition)*Ideology_scaled, data = S3_D3)
S3_reg7_ideo <- lm(mean_National_Pride_scaled ~ factor(Condition)*Ideology_scaled, data = S3_D3)
summary(S3_reg4_ideo)
summary(S3_reg5_ideo)
summary(S3_reg6_ideo)
summary(S3_reg7_ideo)


##############
## Mediation through status threat
##############
library(lavaan)

## Dummy-code conditions
table(S3_D3$Condition)
S3_D3$Condition_CIIM <- ifelse (S3_D3$Condition =="Common Identity", 1, 0)

## Comparison CIIM vs. control
## insert DV as needed
S3_D3_simple_mediation <- '
## Direct effect condition on outcomes (path c) 
Reps_scaled ~ c1*Condition_CIIM

## Indirect effect 1: condition on reported American identity (path a)
mean_Status_Threat_scaled ~ a1*Condition_CIIM

## Indirect effect 2: reported American identity on status threat + controls (path b)
Reps_scaled ~ b1*mean_Status_Threat_scaled

## direct effect
direct_CIIM := c1

## indirect effect (a*b): Sobel test (Delta method)
indirect_CIIM := a1*b1

## total effect
total_CIIM := c1 + (a1*b1)
'

## comparing CIIM vs. control
S3_D3_simple_mediation_Outgroup_fit <- sem(S3_D3_simple_mediation, se = "bootstrap", bootstrap = 5000, data = S3_D3)
S3_D3_simple_mediation_Inclusiveness_fit <- sem(S3_D3_simple_mediation, se = "bootstrap", bootstrap = 5000, data = S3_D3)
S3_D3_simple_mediation_Americanness_fit <- sem(S3_D3_simple_mediation, se = "bootstrap", bootstrap = 5000, data = S3_D3)
S3_D3_simple_mediation_Democrats_fit <- sem(S3_D3_simple_mediation, se = "bootstrap", bootstrap = 5000, data = S3_D3)
S3_D3_simple_mediation_Republicans_fit <- sem(S3_D3_simple_mediation, se = "bootstrap", bootstrap = 5000, data = S3_D3)
summary(S3_D3_simple_mediation_Outgroup_fit, fit.measures = TRUE, rsquare = TRUE)
summary(S3_D3_simple_mediation_Inclusiveness_fit, fit.measures = TRUE, rsquare = TRUE)
summary(S3_D3_simple_mediation_Americanness_fit, fit.measures = TRUE, rsquare = TRUE)
summary(S3_D3_simple_mediation_Democrats_fit, fit.measures = TRUE, rsquare = TRUE)
summary(S3_D3_simple_mediation_Republicans_fit, fit.measures = TRUE, rsquare = TRUE)


##############
## Mediation through status threat and national identity + controls
##############
library(lavaan)

## Dummy-code conditions
table(S3_D3$Condition)
S3_D3$Condition_CIIM <- ifelse (S3_D3$Condition =="Common Identity", 1, 0)

## Comparison CIIM vs. control
## insert DV as needed
S3_D3_serial_mediation <- '
Reps_scaled ~ c1*Condition_CIIM

## Indirect effect 1: condition on reported American identity and controls (path a)
mean_American_Identity_scaled ~ a1*Condition_CIIM
mean_White_Identity_scaled ~ a2*Condition_CIIM
mean_Nationalism_scaled ~ a3*Condition_CIIM
mean_National_Pride_scaled ~ a4*Condition_CIIM

## Indirect effect 2: reported American identity on status threat + controls (path b)
mean_Status_Threat_scaled ~ b1*mean_American_Identity_scaled + b2*mean_White_Identity_scaled + b3*mean_Nationalism_scaled + b4*mean_National_Pride_scaled

## Indirect effect 3: status threat on DV (path d)
Reps_scaled ~ d1*mean_Status_Threat_scaled

## direct effect
direct_CIIM := c1

## indirect effect (a*b): Sobel test (Delta method)
Serial_CIIM_AM_ID := a1*(b1)*d1
Serial_CIIM_WH_ID := a2*(b2)*d1
Serial_CIIM_NA := a3*(b3)*d1
Serial_CIIM_NP := a4*(b4)*d1

## total effect
total_CIIM := c1 + (a1*(b1+b2+b3+b4)*d1) + (a2*(b1+b2+b3+b4)*d1) + (a3*(b1+b2+b3+b4)*d1) + (a4*(b1+b2+b3+b4)*d1)

## covariances
mean_American_Identity_scaled ~~ mean_White_Identity_scaled
mean_American_Identity_scaled ~~ mean_Nationalism_scaled
mean_American_Identity_scaled ~~ mean_National_Pride_scaled
mean_White_Identity_scaled ~~ mean_Nationalism_scaled
mean_White_Identity_scaled ~~ mean_National_Pride_scaled
mean_Nationalism_scaled ~~ mean_National_Pride_scaled
'

## comparing CIIM vs. control
S3_D3_serial_mediation_Outgroup_fit <- sem(S3_D3_serial_mediation, se = "bootstrap", bootstrap = 5000, data = S3_D3)
S3_D3_serial_mediation_Inclusiveness_fit <- sem(S3_D3_serial_mediation, se = "bootstrap", bootstrap = 5000, data = S3_D3)
S3_D3_serial_mediation_Americanness_fit <- sem(S3_D3_serial_mediation, se = "bootstrap", bootstrap = 5000, data = S3_D3)
S3_D3_serial_mediation_Democrats_fit <- sem(S3_D3_serial_mediation, se = "bootstrap", bootstrap = 5000, data = S3_D3)
S3_D3_serial_mediation_Republicans_fit <- sem(S3_D3_serial_mediation, se = "bootstrap", bootstrap = 5000, data = S3_D3)
summary(S3_D3_serial_mediation_Outgroup_fit, fit.measures = TRUE, rsquare = TRUE)
summary(S3_D3_serial_mediation_Inclusiveness_fit, fit.measures = TRUE, rsquare = TRUE)
summary(S3_D3_serial_mediation_Americanness_fit, fit.measures = TRUE, rsquare = TRUE)
summary(S3_D3_serial_mediation_Democrats_fit, fit.measures = TRUE, rsquare = TRUE)
summary(S3_D3_serial_mediation_Republicans_fit, fit.measures = TRUE, rsquare = TRUE)

##############
## Exploratory: Effect on affective polarization
##############
## Compute AP-score
S3_D3$AP_score <- abs((S3_D3$Thermometers_3 - S3_D3$Thermometers_4))
head(S3_D3[1:6,c("Thermometers_3","Thermometers_4","AP_score")])

## Run regressions
S3_aovAP <- aov(AP_score ~ factor(Condition), data = S3_D3)
summary(S3_aovAP)
S3_regAP <- lm(AP_score ~ factor(Condition), data = S3_D3)
summary(S3_regAP)
cohens_f(S3_regAP, partial = FALSE)






##############
## End of the syntax.




